Average TPM Across Time of Select Genes (GSE137001)

# Download data
url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE137001&format=file&file=GSE137001%5FTPM%5Freprogramming%5Fintermediates%2Exlsx"
filename <- "GSE137001_TPM_reprogramming_intermediates.xlsx"
if (!file.exists(filename)) {
  download.file(url, filename)
}
tpms <- read_excel(filename)

tpms_mat <- tpms %>% select(3:41) %>% as.matrix()

metadata <- data.frame(sample_name = colnames(tpms_mat)) %>%
  mutate(group = gsub("(d\\d+_|_R\\d)", "", sample_name)) %>%
  mutate(days = case_when(
    startsWith(sample_name, "d") ~ sub("d(\\d+).+", "\\1", sample_name),
    TRUE ~ "NA"
  )) %>%
  mutate(rep = gsub(".+R(\\d)", "\\1", sample_name)) %>%
  column_to_rownames("sample_name")

PCA

# I don't think this transformation is correct
trans_mat <- log10(1 + (tpms_mat %>% t()))

pca_res <- stats::prcomp(trans_mat)

pca_df <- pca_res$x %>%
  as.data.frame() %>%
  rownames_to_column("sample_name") %>%
  mutate(group = gsub("(d\\d+_|_R\\d)", "", sample_name)) %>%
  mutate(days = case_when(
    startsWith(sample_name, "d") ~ sub("d(\\d+).+", "\\1", sample_name),
    TRUE ~ "NA"
  ))

pca_df %>%
  filter(group %in% c("SKM", "OSKM")) %>%
  ggplot(aes(
    x = PC1, y = PC2, color = days, shape = group
  )) + 
  geom_point(size = 3)

# Get gene list to filter
stress_pathways <- c("GO:0036500")  # Just the ATF6 pathway
# stress_pathways <- c("GO:0009267", "GO:0034198", "GO:0036003", "GO:0036500", "GO:1990928", "GO:0034599", "GO:0034976", "GO:0016236", "GO:0010506")

gostres <- gprofiler2::gost(stress_pathways)
stress_ensg <- gostres$meta$genes_metadata$query$query_1$ensgs

ens2sym <- AnnotationDbi::select(
  EnsDb.Hsapiens.v86,
  keys = keys(EnsDb.Hsapiens.v86),
  columns = c("SYMBOL")
)

stress_genes <- ens2sym %>%
  filter(GENEID %in% stress_ensg) %>%
  pull(SYMBOL)

# Filter data for selected stress genes
tpms$Name <- unlist(lapply(tpms$Name, toupper))  # Capitalize gene names
cat("% stress genes in dataset:", 100*mean(stress_genes %in% tpms$Name))

% stress genes in dataset: 100

stress_tpms <- tpms %>% filter(Name %in% stress_genes)

# Select only the SKM and OSKM samples
# Scale the TPMs of each gene to be between 0 and 1
# Note: I couldn't figure out how to scale row-wise, so I transposed instead
scaled_tpms <- stress_tpms[3:41] %>%
  select(contains("SKM_R")) %>%
  t() %>%
  as.data.frame() %>%
  mutate_all(scales::rescale) %>%  # Min-max scaling to be between 0 and 1
  t() %>%
  as_tibble() %>%
  mutate(Gene = stress_tpms$Name) %>%
  relocate(Gene)

# Format data from wide to long for plotting
df <- scaled_tpms %>%
  pivot_longer(
    contains("_"),
    names_to = c("Day", "Factors", "Rep"),
    names_pattern = "d(\\d+)_(\\w+)_R(\\d+)"
  ) %>%
  rename(scaled_TPM = value) %>%
  mutate_if(is.character, as.factor)

Boxplot

df %>%
  ggplot(aes(x = Day, y = scaled_TPM)) +
  stat_boxplot(geom ='errorbar', width = 0.5) +
  geom_boxplot() +
  facet_wrap(~Factors) + 
  geom_jitter(color = "blue", width = 0.05)

Lineplot

ggp <- df %>%
  ggplot(aes(x = Day, y = scaled_TPM, group = Gene)) +
  geom_line(aes(color = Gene)) +
  facet_wrap(~Factors + Rep)

ggplotly(ggp)

res <- aov(scaled_TPM ~ Day, data = df %>% filter(Factors == "OSKM"))
summary(res)

            Df Sum Sq Mean Sq F value   Pr(>F)    
Day          3  1.946  0.6487   7.727 0.000143 ***
Residuals   76  6.380  0.0840                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res <- aov(scaled_TPM ~ Day, data = df %>% filter(Factors == "SKM"))
summary(res)

            Df Sum Sq Mean Sq F value   Pr(>F)    
Day          3  1.369  0.4562   6.295 0.000715 ***
Residuals   76  5.507  0.0725                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1